
library("AER")

set.seed(123)


simIV <- function(n = 200,
	mu = 0, # exclusion restriction failure
	kappa = 1, # treatment strength
	lambda = 0, # responsiveness
	hetero = 1,
	cov = 0.5,
	sigma = 2 # std for errors
	) {
	# errors
	errors <- MASS::mvrnorm(n, c(1, 1), sigma^2* matrix(c(1, cov, cov, 1), 2, 2))
	coefs <- MASS::mvrnorm(n, c(2, 1), hetero * 0.25 * sigma^2* matrix(c(1, lambda, lambda, 0.5), 2, 2))
	v <- errors[, 1]
	u <- errors[, 2]
	# coefficients
	beta <- coefs[,1]
	pi <- coefs[,2]
	# instrument, treatment, outcome
	z <- rnorm(n, 1, 1)
	delta <- sapply(sapply(kappa * pi, min, 1), max, 0)
	x <- delta * z + (1-delta) * rnorm(n, 0, 1) + 0.2*v
	y <- 5 + beta*x + mu * z +  u + rnorm(n, 0, 2)
	# save
	out <- cbind.data.frame(y, x, z, beta, pi, delta, u, v)
	return(out)
}

runSim <- function(nsims = 10, mu = 0, lambdas, kappas, hetero = hetero, cov = 0.5) {
	results <- matrix(NA, length(lambdas) * length(kappas) * nsims, 7) # save coefs, ratio and rho
	colnames(results) <- c("Coef.OLS","Coef.2SLS","ratio","ratio2","rho","sig.OLS","sig.IV")
	p <- 1
	for (lambda in lambdas) {
		for (kappa in kappas) {
			for (k in 1:nsims) {
				d <- simIV(n = 1000, mu = mu, lambda = lambda, kappa = kappa, hetero = hetero, cov = cov)
				# OLS
				out.ols <- lm(y ~ x, data = d)
				coef.ols <- out.ols$coefficients[2]
				pvalue <- summary(out.ols)$coefficients[2,4]
				sig.ols <- ifelse(pvalue<0.05, 1, 0)
				# IV
				out.iv <- ivreg(y ~ x | z, data = d)
				coef.iv <- out.iv$coefficients[2]
				pvalue <- summary(out.iv)$coefficients[2,4]
				sig.iv <- ifelse(pvalue<0.05, 1, 0)
				# result
				ratio <- abs(coef.iv/coef.ols)
				ratio2 <- abs(coef.iv - coef.ols)/abs(coef.ols)
				rho <- cor(d$x, d$z)
				results[p,] <- c(coef.ols, coef.iv, ratio, ratio2, rho, sig.ols, sig.iv)
				p <- p + 1
				if (p%%100==0) {cat(".")}
			}
		}
	}
	return(results)
}



#########################################################
## No Heterogeneity (Positive Selection)
#########################################################

set.seed(123)
out.c1 <- runSim(nsims = 50, mu = 0, lambda = 0.7, hetero = 0, cov = 0.7,
	kappas = seq(0.1, 1, length.out = 40))

# Exclusion restriction failure
out.c2 <- runSim(nsims = 50, mu = 1, lambda = 0.7, hetero = 0, cov = 0.7,
	kappas = seq(0.1, 1, length.out = 40))

## Heterogeneity (Positive Selection)
out.c3 <- runSim(nsims = 50, mu = 0, lambda = 0.7, hetero = 0.3, cov = 0.7,
	kappas = seq(0.1, 1, length.out = 40))

cat("Done\n")


save(out.c1, out.c2, out.c3, file = "graphs/sim_cont.RData")

#########################################################
## Graphics
#########################################################


load("graphs/sim_cont.RData")
library("latex2exp")
xlab <- TeX(r'($|\hat{\rho}(d, \hat{d})|$)')
ylab <- TeX(r'($|\hat{\tau}_{2SLS} / \hat{\tau}_{OLS}|$)')

# benchmark
coef.ols <- out.c1[,1]
coef.iv <- out.c1[,2]
ratio <- out.c1[, 3]
ratio2 <- out.c1[, 4]
rho <- out.c1[, 5]
sig <- ifelse(out.c1[, 7] == 1, TRUE, FALSE)

pdf("graphs/FigureA6a.pdf")
par(mar = c(4,4,1,1))
plot(1, type = "n", xlim = c(-0.05, 1), ylim = c(-1.5,1.5), axes = FALSE,
		xlab = "", ylab = "")
points(abs(rho), log10(abs(ratio)), col = "#000000", cex = 0.8)
reg <- lm(log10(abs(ratio))~abs(rho))
abline(reg, lwd = 2, col = 2)
y0 <- c(0.001, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100)
axis(1); box()
axis(2, at = log10(y0), labels = y0)
f <- summary(reg)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
text(0.9, 1.45, paste("R2:",sprintf("%.3f",summary(reg)$r.squared)), cex = 1.2)
mtext(xlab, 1, 3)
mtext(ylab, 2, 2.5)
graphics.off()

# publication bias
pdf("graphs/FigureA7a.pdf")
par(mar = c(4,4,1,1))
plot(1, type = "n", xlim = c(-0.05, 1), ylim = c(-1.5,1.5), axes = FALSE,
		xlab = "", ylab = "")
points(abs(rho[sig]), log10(abs(ratio[sig])), col = "#000000", cex = 0.8)
wiv <- which(rho < 0.5 & sig == TRUE)
reg <- lm(log10(abs(ratio[sig])) ~ abs(rho[sig]))
reg2 <- lm(log10(abs(ratio[wiv])) ~ abs(rho[wiv]))
abline(reg, lwd = 2, col = 2) #abline(reg2, lwd = 2, col = 4)
y0 <- c(0.001, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100)
axis(1); box()
axis(2, at = log10(y0), labels = y0)
f <- summary(reg)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
text(0.9, 1.45, paste("R2:",sprintf("%.3f",summary(reg)$r.squared)), cex = 1.2)
mtext(xlab, 1, 3)
mtext(ylab, 2, 2.5)
graphics.off()

# exclusion restriction failure
coef.ols <- out.c2[,1]
coef.iv <- out.c2[,2]
ratio <-out.c2[,3]
ratio2 <-out.c2[,4]
rho <- out.c2[,5]
sig <- ifelse(out.c2[,7]==1, TRUE, FALSE)

pdf("graphs/FigureA6b.pdf")
par(mar = c(4,4,1,1))
plot(1, type = "n", xlim = c(-0.05, 1), ylim = c(-1.5,1.5), axes = FALSE,
		xlab = "", ylab = "")
points(abs(rho), log10(abs(ratio)), col = "#000000", cex = 0.8)
reg <- lm(log10(abs(ratio))~abs(rho))
abline(reg, lwd = 2, col = 2)
y0 <- c(0.001, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100)
axis(1); box()
axis(2, at = log10(y0), labels = y0)
f <- summary(reg)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
text(0.9, 1.45, paste("R2:",sprintf("%.3f",summary(reg)$r.squared)), cex = 1.2)
mtext(xlab, 1, 3)
mtext(ylab, 2, 2.5)
graphics.off()

# HTE
coef.ols <- out.c3[,1]
coef.iv <- out.c3[,2]
ratio <-out.c3[,3]
ratio2 <-out.c3[,4]
rho <- out.c3[,5]
sig <- ifelse(out.c3[,7]==1, TRUE, FALSE)

pdf("graphs/FigureA8a.pdf")
par(mar = c(4,4,1,1))
plot(1, type = "n", xlim = c(-0.05, 1), ylim = c(-1.5,1.5), axes = FALSE,
		xlab = "", ylab = "")
points(abs(rho), log10(abs(ratio)), col = "#000000", cex = 0.8)
reg <- lm(log10(abs(ratio))~abs(rho))
abline(reg, lwd = 2, col = 2)
y0 <- c(0.001, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100)
axis(1); box()
axis(2, at = log10(y0), labels = y0)
f <- summary(reg)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
text(0.9, 1.45, paste("R2:",sprintf("%.3f",summary(reg)$r.squared)), cex = 1.2)
mtext(xlab, 1, 3)
mtext(ylab, 2, 2.5)
graphics.off()

pdf("graphs/FigureA8b.pdf")
subset <- which(sig == TRUE & sign(out.c3[,1])==sign(out.c3[,2]))
par(mar = c(4,4,1,1))
plot(1, type = "n", xlim = c(-0.05, 1), ylim = c(-1.5,1.5), axes = FALSE,
		xlab = "", ylab = "")
points(abs(rho[subset]), log10(abs(ratio[subset])), col = "#000000", cex = 0.8)
reg <- lm(log10(abs(ratio[subset]))~abs(rho[subset]))
abline(reg, lwd = 2, col = 2)
y0 <- c(0.001, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100)
axis(1); box()
axis(2, at = log10(y0), labels = y0)
f <- summary(reg)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=FALSE)
text(0.9, 1.45, paste("R2:",sprintf("%.3f",summary(reg)$r.squared)), cex = 1.2)
mtext(xlab, 1, 3)
mtext(ylab, 2, 2.5)
graphics.off()


